Adjoint Formulation for an Embedded-Boundary 

Cartesian Method 

Marian Nemec,* Michael J. Aftosmis,* Scott M. Mur man.* 
and Thomas H. Pulliam^ 

NASA Ames Research Center 
MS T27B . Moffett Field . CA 94035 


Abstract submitted to the 43rd AIAA Aerospace Sciences Meeting 
Jan. 10-13, 2005/Reno, NV 

I. Introduction 

M ANY problems in aerodynamic design can be characterized by smooth and convex objective functions. 

This motivates the use of gradient-based algorithms, particularly for problems with a large number 
of design variables, to efficiently determine optimal shapes and configurations that maximize aerodynamic 
performance. Accurate and efficient computation of the gradient, however, remains a challenging task. In 
optimization problems where the number of design variables dominates the number of objectives and flow- 
dependent constraints, the cost of gradient computations can be significantly reduced by the use of the 
adjoint method. 

The problem of aerodynamic optimization using the adjoint method has been analyzed and validated 
for both structured and unstructured grids. The method has been applied to design problems governed 
by the potential, Euler, and Navier-Stokes equations and can be subdivided into the continuous 1,2,3,4 and 
discrete formulations . 0 ’ 6 ’ ' ,8j9, 10,11 Giles and Pierce 12 provide a detailed review of both approaches. Most 
implementations rely on grid-perturbation or mapping procedures during the gradient computation that 
explicitly couple changes in the surface shape to the volume grid. The solution of the adjoint equation is 
usually accomplished using the same scheme that solves the governing flow equations. Examples of such 
code reuse include multistage Runge-Kutta schemes coupled with multigrid , 3 ’ 9 approximate-factorization , 6 
line-implicit Gauss-Seidel , 8, 13 and also preconditioned GMRES . 10,11 

The development of the adjoint method for aerodynamic optimization problems on Cartesian grids has 
been limited. In contrast to implementations on structured and unstructured grids, Cartesian grid methods 
decouple the surface discretization from the volume grid . 14 This feature makes Cartesian methods welL 
suited for the automated analysis of complex geometry problems, and consequently a promising approach 
to aerodynamic optimization. Melvin et al 15 developed an adjoint formulation for the TRANAIR code , 16 
which is based on the full-potential equation with viscous corrections. More recently, Dadone and Gross- 
man 17 presented an adjoint formulation for the Euler equations. In both approaches, a boundary condition 
is introduced to approximate the effects of the evolving surface shape that results in accurate gradient 
computation. 
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In Ref. 18, we presented an effective framework for complex geometry aerodynamic optimization problems. 
Key modules of the framework were a. Cartesian inviscid-flow anafysis package, CartSD . 19:20 and a direct- 
CAD interface. 21 A genetic and gradient-based algorithms were used to drive the optimization procedure. 
The gradients were computed via finite-differences. 

In this work, we improve the efficiency of the gradient- computation by developing an adjoint method for 
aerodynamic optimization on Cartesian grids with embedded boundaries. We present the numerical imple- 
mentation for the discrete approach, where the adjoint solver efficiently reuses the domain decomposition, 
multigrid, and time-marching schemes of the flow solver. We focus on the accuracy of the gradient computa- 
tion, which requires careful treatment of cut-cells, the underlying surface triangulation, and the interface to 
a parametric-CAD system. The final paper will contain a number of complex geometry, industrially relevant, 
examples with many design variables to demonstrate the effectiveness of the adjoint method on Cartesian 
grids . 

II. Problem Formulation 

The aerodynamic shape optimization problem consists of 
determining values of design variables X , such that the ob- 
jective function J is minimized 

min J(X,Q) (1) 

subject to constraint equations Cj 

Cj(X, Q) <0 j = l ,...,N e (2) 

where the vector Q denotes the continuous, conservative flow 
variables and N c denotes the number of constraint equations. 

The flow variables are forced to satisfy the governing flow' 
equations within a feasible region of the design space Q 

T{X,Q) = 0 Vlefl (3) 

which implicitly defines Q — f(X). 

The objective function defines the goals of the optimiza- 
tion problem, while the constraint equations limit the feasi- 
ble region of the design space. The constraints may involve 
performance functionals, such as lift, geometric quantities, 
such as volumes and thicknesses, and also simple bound con- 
straints for design variables. Performance objectives can be 
specified for the entire configuration or for a specific sub- 
set of components, for example moments on control surfaces. ^ , . , , 

r c _ Figure 1. ComDonents of the analysis module 

This is accomplished by using the Geometry Manipulation 

Protocol, 22 where we specify component hierarchies for parametric-CAD assemblies and triangulations to 
intuitively reflect aerodynamic design goals. 

A modular framework 18 is used to solve the optimization problem defined by Eqs. 1-3. We cast the 
optimization problem as an unconstrained problem by lifting the side constraints, Eq. 2, into the objective 
function using a penalty method. The constraint imposed by the flowfield equations, Eq. 3, is satisfied at 
every point within the feasible design space, and consequently these equations do not explicitly appear in the 
formulation of the optimization problem. An unconstrained BFGS quasi-Newton algorithm coupled with a 
backtracking line search 23,24 is used to find the optimal solution. At the core of the optimization framework 
is the analysis module, w-hich consists of a CAD system interface provided by CAPRI. 25,21 a Cartesian grid 
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generator for component-based geometry, 19 and a flow solver, as outlined in Fig. 1. Below, we provide a 
brief description of the flow solver. Thereafter, we focus on the development of the gradient computation 
algorithm. 


III. Flow Solver 

The governing flow equations are the three-dimensional steady-state Euler equations of a perfect gas. A 
second-order accurate finite-volume spatial discretization is used, which is based on van Leer’s flux vector 
splitting. 20 The resulting system of nonlinear equations is given by 

R(Q,X) = o (4) 

where Q ~ [p, pu, pv, pu, pE] T denotes the discrete vector of integral cell averages. In each cell of the volume 
grid, the residual vector can be expressed as 

^ 41 ^“ ( 5 ) 

u da 

where Vq is the cell volume, dQ. are the faces of the cell, F is the flux vector, and n is the face unit normal. 
Although R is written as a function of the design variables, we emphasize that during a flow solution the 
design variables remain constant. A solution is achieved by introducing an unsteady term given by 

^ + R(Q) = 0 (6) 

and marching to steady-state using a Runge-Kutta scheme in conjunction with a parallel multigrid method. 20 


IV. Adjoints and Sensitivities 


The gradient, G, of the objective function J [X, Q(X)] is given by 

dj ^^7 , dJdQ 
T dX 8X ^ dQ dX 


(7) 


where Q denotes a discrete steady-state solution of the flow equations. We reduce the vector of design 
variables, X , to a scalar in order to clearly distinguish between partial and total derivatives. For problems 
with multiple design variables, it may be helpful to note that G and dJ/dX are [1 x JVp] row vectors, 
dJ/dQ is a [1 x JV>] row vector, and dQ/dX is a [iV F x iV D ] matrix, where No and iYp represent the number 
of design and flow variables, respectively. 

Throughout this development, we assume that the implicit function Q(X) is sufficiently smooth. The- 
oretically, this assumption is violated at flow- discontinuities, for example shock waves, and leads to the 
requirement for additional boundary conditions. In practice, the smoothing of the discontinuities by the 
underlying spatiai discretization of the flow equations provides a Q^X'j that is sufficiently smooth. 

In Eq. 7, the evaluation of the term dQ/dX, referred to as the flow sensitivities, is obtained by differen- 
tiating Eq. 4 with respect to the design variables 


dR _ OR ( dRdQ 
dX ~ dX + dQ dX 


( 8 ) 


Realizing that ^ = 0, since for any design variable Eq. 4 is always satisfied, Eq. 8 simplifies to the following 
large sparse system of linear equations 

dR dQ _ dR 

dQ dX ~ dX 


(9) 
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The direct, or flow-sensitivity, method results from solving Eq. 9 for the flow sensitivities dQ/dX and using 
these values in Eq. 7 to obtain the gradient. 

In order to formulate the discrete- adjoint method, substitute Eq. 9 into Eq. 7 to obtain 


dj _ (dR y 1 dR 

dX ~ dX dQ\dQj dX 

From the triple-product term in Eq. 10, define the following intermediate problem 

, T w f^Y 1 

dQ { dQj 


( 10 ) 


(ii) 


w'here ip is a [iV> x 1] column vector. Post-multiplication of both sides by dR/dQ and applying the transpose 
operator results in the following linear system of equations 


dY , dY 

dQ * ~ dQ 


( 12 ) 


This is known as the adjoint equation, and the vector ip represents the adjoint variables. Substituting ip into 
Eq. 10, the expression for the gradient becomes 


d J __ dj , T dR 

dX = dX ~ v dX 


(13) 


Note that Eq. 9 is a linear system with multiple right-hand sides dependent on the number of design variables. 
In contrast, Eq, 12 is independent of the design variables, which is the reason for the efficiency of the adjoint 
method. 


A. Numerical Method 

The parallel multigrid method from the flow solver is adopted to solve both the adjoint and flow-sensitivity 
equations. We consider the flow sensitivity equation first, since with this approach the explicit linearization 
of the residual equations is not necessary. The matrix-vector product on the left-hand side of Eq. 9 is 
approximated as follows. 

dR R(Q + tv)-R(Q) 

9Q 7 (14) 

wffiere v — dQ/dX and e = \/( € m)/IMI- F° r cases considered in this work, we set e = 1 x 10~ 4 in order 
to simplify the numerical implementation of Eq. 14. To solve for the flow-sensitivities, let 

Q* = Q -f ev (15) 

and introduce an unsteady term in Eq. 9 to obtain the following expression 

^ + - «(«■)=■ *(«)-'§ (i6 > 

The only difference between Eqs. 16 and 6 is an extra right-hand-side source term. Recall that Q is held 
fixed and we set the initial vector v to zero. Hence, the solution algorithm used for the flow' equations can 
be recycled to solve the flow-sensitivity equation with only a minor modification to the multigrid forcing 
function. We expect the same asymptotic convergence rate for both algorithms. A limitation of this approach 
is the choice of the constant stepsize e. In our experience, the approach works well for subsonic flows and 
serves as a useful debugging tool for establishing the accuracy of the gradient computation. The final paper 
will include a detailed discussion of the adjoint solver. 
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The partial derivative term dJ/dQ in Eqs. 7 and 12 is evaluated using finite differences. Note that care 
must be taken to ensure that the evaluation of the objective function is consistent with the enforcement of 
the boundary conditions in the residual equations. Otherwise, the adjoint solution may contain oscillations 
near wall boundaries, as demonstrated by Lu and Darmofal. 2, 'We plan to investigate this further in the 
final paper. 

The remaining partial derivative terms in Eqs. 7 and 13, namely the objective function sensitivity dJ/dX 
and the residual sensitivity dR/dX , are usually approximated with finite differences as well. Although in 
some instances the partial derivatives can be obtained analytically, finite differences are necessary when a 
CAD system is used to parameterize the geometry'. For example, the centered- difference formula for the 
residual sensitivities is given by 


dR RjX + he n , Q)-R (X - he n , Q) 

dX n ~ 2/i ( () 

where e n denotes the nth unit vector, 

h = max (e 0 • |X n |, v^m) ( 18 ) 

and n = 1, . . . , Aq. Typical value of e 0 is lx 10“ 3 . It is important to realize that Eq. 17 involves two 
evaluations of only the residual vector per design variable and not two flow solutions, i.e. Q is constant. 

Before presenting our numerical implementation for the computation of dJ/dX and dR/dX on Cartesian 
grids, it is insightful to briefly review the implementations on structured and unstructured grids. The 
standard finite-difference approach relies on a grid-perturbation strategy to smoothly “convect” the nodes 
of the volume grid with the perturbed surface geometry. This results in an accurate and efficient evaluation 
of the partial derivative terms. Alonso et al 28 have successfully extended this approach to a CAD-based 
design environment by introducing geometry patches within the parametric-CAD model. 

An important observation is that the grid-perturbation strategy explicitly couples shape sensitivities to 
volume grid sensitivities. However, a grid-perturbation without a corresponding surface perturbation should 
not, in principle, influence the objective function gradient. On the basis of this argument, Jameson and 
Kim 29 introduced a reduced adjoint gradient formulation where only surface perturbations are considered. 
Their results indicate that for shape optimization problems the reduced approach works well. Similar results 
are also presented by Soto and Lohner. 30 In contrast, work by Anderson and Venkatakrishnan 4 and later 
by Nielsen and Anderson 31 shows that for cases that involve rigid body' motion, for example flap position 
optimization, the use of a grid-perturbation strategy is important. This appears to be related to the treatment 
of the trailing-edge singularity. 

Turning our attention back to embedded-boundary Cartesian grids, we note that a perturbation of the 
surface shape affects only r a few* near-by cells of the volume grid. Therefore, the extent of grid sensitivities 
is naturally limited to roughly 0(N 2 ) cells for a volume grid with 0('N 3 ) cells. However, a straightforward 
implementation of finite differences is complicated by the emergence and disappearance of cut-cells at the 
w'all boundary, as well as the relative motion of the discretized surface with respect to the volume grid. 

To circumvent these difficulties, we use the following approach. For dj jdX, the converged flow solution 
is reconstructed to the vertices of the underlying surface tri angulation. Given a small perturbation, a 
new surface triangulation is generated using the CAPRI interface. We force the new triangulation to have 
a one-to-one triangle mapping with respect to the baseline triangulation. This is accomplished by using 
the parametric values of the baseline triangulation wffien tessellating the perturbed surface. The objective 
function values required for the finite-difference approximation are now evaluated using the triangulations 
instead of the volume grid. Hence, the potentially non-smooth changes in the cut-cells are avoided, and 
the enforcement of constant Q in the evaluation of this partial derivative term becomes a trivial task. This 
approach works well for both shape and rigid-body motion design variables. 

For the remaining term dR/dX, we have implemented and studied two approaches, both based on finite- 
differences. In the first approach, we assume that the only parameter dependent on the design variables 
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in Eq. 4 is the surface normal. This approach is similar to the "transpiration" boundary condition used in 
TRAN AIR.. 32 - 15 It is relatively easy to implement, and the evaluation of 8R/dX is fast since the volume grid 
remains unchanged. The accuracy of this approximation can be sufficient when the design variables involve 
shape changes normal to the surface including wing-twist distributions. However, a possible limitation of 
this approach may be rigid-body motion. 

In the second approach, w-e construct a new r volume grid for the perturbed surface triangulation. The 
baseline solution is interpolated to the new' grid via a fast interpolation algorithm presented in Ref. 33. We 
evaluate the residual on the new' grid, w’hich is then passed back to the baseline grid for the finite-difference 
computation. Cut-cells that appear or disappear during this procedure are tagged, and we either zero their 
contribution to the residual sensitivity computation, or we interpolate a value from their face neighbors. We 
evaluate both approaches in the following section. 

V. Results and Discussion 

In this abstract we study a simple, two-dimensional design example to validate our approach. The baseline 
geometry is the NACA 0012 airfoil parameterized with a cubic B-spline curve using 17 control points. We 
choose the vertical motion of a control point near the leading edge of the airfoil as the design variable, shown 
in Fig. 2. The freestream Mach number is 0.5 and the angle of incidence is fixed at 2 deg. The goal of the 
optimization is lift enhancement. The baseline Cl is 0.245 and we specify a target Cl of 0.5. The “volume” 
grid contains roughly 9, 000 cells, but since this is a two-dimensional problem, the depth in the span- wise 
direction contains only two cells. To help validate the results, w^e solve a similar problem using Optima2D, n 
which is an aerodynamic optimization tool based on a structured-grid approach. 

We begin by examining the computation of the 
dj jdX term in Eq. 7. We evaluate the partial 
^grivative of the lift fun c t i o n al (objective func- 
tion) for the above problem with the angle of inci- 
dence set to 0 and 2 deg. Note that at 0 deg., this 
partial derivative term should vanish, providing a 
good validation exercise. The result at 0 deg. is 
indeed zero, and at 2 deg. we obtain —2.81 x 10 -3 . 

Optima2D predicts — 1.4 x 10~ 2 , which is a rea- 
sonable agreement considering the differences of 
the tw'o approaches. 

Prior to solving the flow-sensitivity equa- 
tion, we evaluate the residual sensitivity dR/d X. _ . . , . . _ ^ 

* __ Figure 2. Design variable and B-spline parameterization 

Figure 3(a) shows dR/d A values for the y- of the NACA 0012 airfoil 

momentum equation. These are obtained by con- 
structing (or re-cutting) a new volume grid for the perturbed shape and interpolating the steady-state 
solution. Figure 3(b) shows the results wffien only the surface normals are changed. The agreement in Fig. 3 
is excellent; however, we note that for the remaining field variables significant errors may arise due to this 
approximation. W r e plan to investigate this further in the final paper. 

The solution of the flow’-sensitivity equation is shown in Fig. 4. Note that the zone of influence for 
the shape design variable extends from the leading edge to roughly 60% chord on the upper surface. The 
convergence of the flow and sensitivity equations is shown in Fig. 5. The flow equations converge in 500 
multigrid cycles (only a 2-level multigrid is used), and the flow sensitivities require additional 400 cycles. The 
final gradient value is —1.35. Optima2D computes a value of —1.67, while a finite- difference estimation of 
the gradient by recomputing the flow gives —1.55. A grid-refinement, study will be included in the final paper 
to demonstrate gradient accuracy. Furthermore, the final paper will contain additional design examples to 
demonstrate the effectiveness of the proposed approach. 
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Figure 3. Contours of dR/dX for the y-momentum equation 
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(a) Cartesian- grid formulation 


(b) Structured-grid approach using Optimal D 11 

Figure 4. Contour plots of flow sensitivity: changes in density with respect to a shape design variable 
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Figure 5. Convergence to steady-state for the mean flow and the flow-sensitivity equation 
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